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SECTION  I 


INTRODUCTION 

Even  though  finite  difference  fluid  mechanics  techniques  have 
matured  significantly  in  recent  years,  the  "computational  wind  tunnel" 
remains  a  distant  vision.  Physically  complete  flow  field  predictions 
have  been  generally  limited,  by  computer  speed  and  cost  (at  least),  to 
relatively  small  regions  on  two-dimensional  space.  However,  with  the 
current  rapid  evolution  of  digital  hardware  technology,  the  near-term 
prospect  exists  for  emergence  of  a  large-scale  computational  capability 
for  aerodynamic  flowfield  prediction.  Candidate  fourth  generation  com¬ 
puters  include  at  least  the  STAR,  and  CYBER,  and  the  CRAY,  and  perhaps 
a  dedicated  Navier-Stokes  solver  (Proceedings  of  the  NASA  Workshop 
on  Future  Computer  Requirements  for  Computational  Aerodynamics,  ref.l). 

Numerical  algorithms  for  solution  of  simplified  forms  of  the  three- 
dimensional  Navier-Stokes  equations  have  historically  been  based  upon 
finite  difference  theory,  (ref.  2).  MacCormack's  explicit 
time-split  method,  MacCormack  and  Paullay  (ref.  3),  emerged  as  an 
efficient  second-order  accurate  algorithm  for  solution  of  the  inviscid 
Euler  equations,  and  combined  with  shock-capturing  has  enjoyed  world-wide 
use.  Alternatively,  two-dimensional  viscous  flow  computations  have 
generally  employed  alternating-direction,  implicit,  second-order  accurate 
finite  difference  algorithms  for  time  (or  space)  evolution,  and  suc¬ 
cessive  overrelaxation  for  boundary  value  problems  and/or  direct  iteration 
at  steady-state.  The  Crank-Nicolson  algorithm  is  the  representative  exam¬ 
ple  procedure  (ref.  2).  An  implicit  algorithm  circumvents  the  severe 
stability  constraint  associated  with  explicit  integration  procedures;  for 
example.  Beam  and  Warming  have  developed  implicit  forms  for  the  Euler  equa¬ 
tions  (ref.  4)  and  the  compressible  Navier-Stokes  equations  (ref.  5). 

In  recent  years,  the  engineering  community  has  become  cognizant  of 
functional  analysis  and  interpolation  theory  as  an  alternative  theoretical 
foundation  for  development  of  numerical  algorithms  for  the  Navier-Stokes 
equations.  The  finite  element  method  is  the  engineering  embodiment  of 
these  tools  of  the  mathematician.  The  apparent  fundamental  difference 
between  finite  difference  and  finite  elemert  theory  is  that  the  former  deals 
with  difference  algebra  while  the  latter  employs  vector  field  theory  and 
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calculus.  However,  as  must  occur,  both  concepts  belong  to  interpolation 
theory,  Prenter  (ref.  6).  A  strong  common  theoretical  foundation  is  thus 
emerging,  which  is  recognized  and  employed  herein  to  assist  in  derivation 
of  optimal ly-accurate,  and  efficient  algorithms  for  the  highly  non-linear 
Navier-Stokes  equations.  For  example,  use  of  a  linear  finite  element  al¬ 
gorithm  applied  to  the  diffusion  term  in  the  Navier-Stokes  equations  on  a 
uniform  grid  yields  identically  the  second-order  accurate  finite  difference 
form,  i.e.,  the  central  difference  operator.  This  does  not  necessarily 
extend  to  initial-value  and/or  non-linear  differential  operators  however. 
For  example,  finite  element  theory  renders  the  linear  finite  element  dis¬ 
crete  embodiment  of  the  substantial  time-derivative. 


a  fourth-order  accurate  representation  on  a  uniform  grid.  Figure  la-b 
illustrates  the  accuracy  obtainable  for  implicit  solution  of  the  (hyper¬ 
bolic)  continuity  equation  using  linear  finite  element  theory.  An  elemen¬ 
tary  rearrangement  produces  the  Crank- Ni col  son  finite  difference  form, 
which  is  only  a  second-order  accurate  algorithm.  The  computed  results  shown 
in  Fig.  lc  illustrate  the  significant  loss  in  accuracy  attendent  with  this 
algorithm  for  this  discretization. 

This  elementary  comparison  simply  documents  that  the  finite  element 
embodiment  of  linear  interpolation  theory  can  produce  optimal ly-accurate 
finite  difference  equivalents  of  differential  operators,  representative  of 
those  dominating  the  three-dimensional  Navier-Stokes  equations.  Of  perhaps 
greater  importance,  finite  element  theory  also  yields  a  thoroughly  standard¬ 
ized  procedure  for  generation  of  potentially  higher-order  accurate  operators 
using  quadratic,  cubic  and/or  higher-dimensional  interpolation  theory. 
Confirmation  for  linear  hyperbolic  and  parabolic  equations  is  presented  in 
reference  7. 

The  current  objective,  for  which  results  of  the  Phase  I  effort  are  reported 
herein,  is  to  assess  extension  of  this  highly  promising  theoretical  procedure 
to  solution  of  the  highly  non-linear,  three-dimensional  Navier-Stokes  equa¬ 
tions  for  aerodynamic  configurations.  The  primary  requirement  of  Phase  I 
was  to  assess  factors  controlling  accuracy,  convergence  and  efficiency  of  the 
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predictions  produced  by  the  developed  finite  element  theory  for  partial 
differential  equations  modeling  the  important  non-linear  convection 
operator  in  the  Navier-Stokes  equations.  The  specific  approach  has  been 
to  delete  all  physical  viscosity  terms,  which  yields  the  Euler  simplifi¬ 
cation  to  Navier-Stokes,  to  determine  algorithm  performance  in  the  ab¬ 
sence  of  the  smoothing  introduced  by  viscosity.  This  corresponds  as  well 
to  the  large  Reynolds  number  limit  of  Navier-Stokes,  wherein  the  field  of 
aerodynamics  generally  resides.  Both  linear  and  quadratic  finite  element 
tensor  product  formulation  have  been  evaluated  and  proven  viable,  as 
presented  herein. 
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SECTION  II 


THEORETICAL  ANALYSIS 


1.  Overview 

The  multi-year  goal  is  to  derive,  implement  into  a  computer 
program,  and  evaluate  accurate  and  economical  numerical  algorithms  for 
solution  of  the  three-dimensional  Navier-Stokes  equations  for  turbulent 
aerodynamic  flows.  This  section  introduces  the  topic  with  a  concise 
statement  of  the  appropriate  differential  equation  system  requiring  solu¬ 
tion.  The  uniformity  pervading  this  system  then  yields  a  concise  state¬ 
ment  of  the  dissipative-Weighted  Residuals  finite  element  algorithm.  A 
theoretical  analysis  of  a  linearized  model  equation  is  discussed  that 
isolates  and  estimates  acceptable  levels  for  key  parameters  imbedded 
within  the  basic  algorithm  that  control  accuracy  and  convergence. 


2.  Navier-Stokes  Equation  System 


The  set  of  partial  differential  equations  governing  the  transient 
three-dimensional  flows  of  interest  is  the  familiar  and  very  non-linear 
Navier-Stokes  system.  All  equation  systems  studied  are  derived  directly 
from  the  Navier-Stokes  equations.  The  non-dimensional  conservation  form 
(in  Cartesian  tensor  notation  with  summation  implied  over  repeated  sub¬ 
scripts)  for  a  compressible,  viscous,  heat-conducting  fluid  is 


L(p)  “ft  +  jljM*0 

3(puJ  a  r  1 

=  -5t-*srb  p“i  •Td'° 

J 

L(pH)  •  *  377[uj  pH  •  “j(TU  +  pSlj>  '  qj]  •  H  ■  0 

J 


(1) 

(2) 

(3) 


In  equations  1  -  3,  p  is  density,  pu..  is  the  momentum  vector,  p  is 
pressure,  and  the  Stokes  stress  tensor  x. .,  heat  flux  vector  qi9  and  stag 

1 J  J 

nation  enthalpy  H  are  defined  as 
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Herein,  y  is  the  dynamic  viscosity,  k  is  thermal  conductivity,  T  is  tem¬ 
perature,  h  is  static  enthalpy,  6..  is  the  Kronecker  delta,  and  Re  is  the 

'  J 

Reynolds  number.  An  equation  of  state  closes  the  system. 

Equations  1-6  are  valid  for  both  laminar  and  turbulent  flows. 

For  the  latter,  however,  their  solution  becomes  tractable  in  a  practical 
sense  only  after  time-averaging.  While  more  definitive  resolutions  are  the 
subject  of  current  research,  conventional  mass-weighted  time-averaging  will 
probably  serve  the  aerodynamic  requirements.  Therefore,  the  Reynolds  de¬ 
composition  yields  (ref.  9) 

u j (  x,t)  =  lTj- (  x)  +  Uj(  x,t)  (7) 

The  mass-weighted  time-average  velocity  is  defined  as 


•  (B). 

(9) 


pu^.  =  pITjU'j  +  PuTuJ  (10) 

Equations  7-8  are  also  employed  to  define  the  time-averaged  and  fluc¬ 
tuating  components  of  the  scalar  fields  h,  H,  and  T.  Hence,  for  example. 
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Substitution  of  the  defining  equations  7-8  into  1-6,  time 
averaging  and  collecting  terms  yields  the  time-averaged  Navier-Stokes 
equation  system 


is  constituted  solely  of  viscous  contribution  to  the  state  of  stress.  The 
Reynolds  stress  tensor  puTuT  represents  six  additional  unknowns  that  have 

*  J 

been  introduced  into  equations  12  -  14.  Various  techniques  are  available 
to  accomplish  closure,  references  8  -  10.  The  primary  emphasis  in  the 
current  assessment  corresponds  to  the  identical  vanishing  of  both  stress 
tensors. 


The  Dissipative-Weighted  Residuals  Algorithm 


The  time-averaged  Navier-Stokes  equations  are  established.  This  system 
encompasses  the  Euler  equations  for  inviscid  flows  as  a  subset,  as  well  as 
the  various  simplified  forms  including  parabolic  Navier-Stokes,  boundary 
layer,  and  incompressible  flows  including  the  streamfunction-vorticity  form. 
Independent  of  the  particular  sub-class,  each  of  the  partial  differential 
equations  of  the  various  coupled  systems  is  of  the  form 


8q. 

L(q3>  ■  ^ 


L  [VJ  +  ^  ♦  M,J  - 


In  equation  16,  q.  is  a  generalized  dependent  variable  (which  could  be  a 
vector,  e.g.,  puj),  u.-  is  the  convection  velocity,  is  a  stress  tensor 
(if  present),  and  f  is  any  non -homogeneity.  The  n-dimensional  problem  is 
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defined  on  the  Euclidean  space  R  ,  spanned  by  the  x  coordinate  system  with 
scalar  components  ,  1  <  i  <_  n,  and  t  is  a  generalized  ini tial -val ue 
coordinate.  The  solution  domain  ft  is  defined  as  the  product  of  Rn  and  x, 
for  all  elements  of  x  belonging  to  Rn  and  all  elements  of  x  belonging  to  the 
open  interval  measured  from  x0,  i.e.. 


ft  =  Rn  x  x  =  {(x,t):  xeRn  and  xe[xQ,x)} 

The  boundary  an  of  the  solution  domain  is  the  product  of  the  boundary  3R  of 
Rn,  spanned  by  x,  and  x,  i.e.,  3ft  h  aR  x  x.  Thereupon,  a  differential  con¬ 
straint  may  be  applied  of  the  form 

f(qj)  =  +  a2  ^  q^f\.  +  a3  =  0  (17) 

In  equation  17,  the  a^  are  specified  coefficients  and  n^  is  the  outwards 
pointing  unit  normal  vector.  Finally,  an  initial  distribution  of  q.  on 
ft0  =  R  x  x0 


q(x,x0)  =  q0(x) 


The  basic  concept  of  a  (finite  element)  numerical  solution  algorithm  is 
to  subdivide  ft  into  non-overlapping  sub-domains  ft  ,  such  that  their  union 
forms  ft,  i.e.. 


ft  =  Ufte  -  URg  x  x  (19) 

and  to  enforce  approximate  satisfaction  of  equations  (16)  -  (18)  on  Ufte>  The 

finite  element  approach  is  to  assume  every  member  of  the  set  q.  interpolated 

on  ft  as 
e 


qQ(x,x)  =  {Mx)}T{Q(t}} 


In  equation  (20),  the  elements  of  the  row  matrix  {Nk(x)}T  are  typically  assumed 
polynomials  on  t,  complete  to  degree  k,  and  the  elements  of  {Q(x)}e  are  the 
(unknown)  expansion  coefficients  that  are  required  determined  by  numerical 
procedures.  The  resultant  solution  q*,  which  is  the  numerical  approximation  to 
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im. 


r 


q(x,r),  is  then  the  summation  of  the  individual  qg  over  the  total  number 
of  elements  (M)  spanning  n  as 

M 

q(x,T)aq*(x,T)  =  I  qp(x,i )  (21) 

e=l  e 

The  fundamental  requirement  in  formulation  of  a  numerical  solution 
algorithm  is  to  render  the  "distance"  between  q*  and  q  a  minimum,  i.e.,  to 
make  the  error  in  the  approximate  solution  q*  as  small  as  possible. 

Since  the  approximate  solution  q*  lies  in  the  domain  spanned  by  the 
interpolation  basis  > ,  then  the  error  can  be  minimized  by  requiring 
it  to  be  orthogonal  to  this  space.  The  measurable  error  is  L(q*)  f  0, 
as  results  from  substitution  of  equation  21  into  16.  From  varia¬ 
tional  boundary  value  theory  (ref.  11),  this  error  is  minimized  (in  an 
appropriate  norm)  by  requiring  it  to  be  orthogonal  to  -CNk > .  Quite  simply, 
then. 


(N.  (x)JL(q*)  =  (0)  (22) 

Rn 

is  the  desired  statement.  (It  is  of  considerable  interest  that,  with  exten¬ 
sion  of  the  "weights"  { Nk )  to  include  other  specified  parameters,  equation 
22  can  be  a  theoretical  basis  for  derivation  of  most  finite  difference  algor¬ 
ithms,  reference  12.  However,  selection  of  weights  other  than  { Nk }  normally 
does  not  enhance  solution  accuracy,  as  discussed  herein,  and  in  reference  7.) 
4n  addition  to  equation  22  is  desired,  as  becomes  apparent  from  the  theoreti¬ 
cal  analysis  in- the  ne*t  section,  since  equation  22  can  yield  a  neutrally 
stable  algorithm  with  no  inherent  dissipative  mechanisms  for  damping  generated 
numerical  error  The  desired  form  can  be  achieved  by  requiring  the  gradient 
of  L(q*)  to  also  be  orthogonal  to  {N^}  ,  i.e., 

6*  jRn{Nk}VL(q*)  =  {0}  (23) 

where  £  is  an  n-dimensional  vector  to  be  determined.  Finally,  for  consis¬ 
tency,  the  error  on  must  also  be  minimized;  hence, 

x  [  {N.H(q*)  =  {0}  (24) 

*  aR 

where  X  is  a  scalar  to  be  determined. 


4 

3 


-9- 


The  dissipative  Weighted-Residuals  algorithmic  statement  is  the 
linear  combination  of  equations  22  -  24.  Since  the  {Nk(x)}  are  defined 
non-zero  only  locally  on  each  subdomain  (finite  element)  £2  ,  the  global 
integrations  indicated  are  replaced  by  integration  on  R^,  and  their  indi¬ 
vidual  contributions  appropriately  summed  ("assembled")  to  yield  the  algo¬ 
rithmic  statement 


(Nk>L(qe)  -  B  •  (Nj(}VL(qe)  +  X 


(NkH(qe) 

3Ra  u  3Rn 
e 


=  (0) 


(25) 


In  equation  25,  Sg  is  the  defined  assembly  operator  which  amounts  to 
row-wise  element  summations  within  the  defined  element  matrices.  Similarly, 
L(q*)  etc.  has  become  replaced  by  an  element  of  its  definition  L(qg),  see 
equation  21. 

Since  the  matrix  elements  of  (N^ }  and  the  various  errors  L(qe),  VL(qe) 
and  £(qg)  all  possess  known  functional  dependence  on  x,  the  indicated  in¬ 
tegrals  are  readily  evaluated.  The  sole  unknown  dependence  is  on  the  ini¬ 
tial  value  coordinate  t;  hence,  equation  25  is  a  system  of  ordinary 
differential  equations  which  express  the  x-derivative  of  q*(x,T)  in  terms 
of  its  current  distribution  of  nodal  values  (Q( t) } .  Therefore,  equation  25 
takes  the  form  of  the  matrix  ordinary  differential  equation  system. 

Se[[C]e(9i;  *  [U]eJ0;e  t  [K]e{Q)e  Mf)J  3.  {0}  .  ..  (26) 

The  terms  in  equation  26  correspond  one-to-one  with  those  in  equation  16, 
i.e.,  the  first  accounts  for  evolution,  the  second  convection,  the  third  the 
stress  field,  and  the  fourth  contains  all  non-homogeneities.  The  [U]  and 
[K]  matrices  in  equation  26  are  general  non-linear  functions  of  { Q0 > ; 
only  the  dominant  non-linearity  is  explicitly  expressed. 

These  theoretical  proceedings  have  thus  transformed  a  single  non-linear 
partial  differential  equation  into  a  much  larger  system  of  non-linear  or¬ 
dinary  differential  equations.  The  requirement  is  to  integrate  this  system. 
For  efficiency,  a  single-step  implicit  algorithm  is  preferred,  and  all  are 
expressed  within  the  expression 
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T - . 

4 


{Q}j+1  5  Wj  +  h[0{Q}j+1  +  0  -  0)tQ>jJ 


(27) 


where  0  <_  8  <_  1  controls  implicitness,  stability  and  accuracy,  and  j  is  the 
integration  step  index  and  h  is  the  step  size.  Equation  26  provides  the 
(implicit)  expression  for  the  {Q}'  column  matrices.  While  not  acceptable 
in  practice,  equations  26  -  27  can  be  combined  using  a  formal  inverse  operation, 
yielding 


(Q}j+1  -  (QJj  -  heECDJi, 


f[u] 


j+i  +  [K]j+i 


{Q}j+i  +  {f}j+i 


h(i  -  eHc]"1 


[Ulj 


tK]d) 


(Q>,  +  if), 


(28) 


Combining  like  terms  and  letting  I  signify  the  identity  matrix,  equation  28 
can  be  rewritten  as 


-  I  -  h(l  -  9)[C]-'((U].  *  [Kjjj 

-  h(e[C]ji,tf}j+1  *  (1  -  SltClj'Wj)  (29) 

which  is  a  non-linear  algebraic  equation  system  for  (Q}j+j* 

Therefore,  the  final  algorithmic  requirement  is  to  establish  a  suitable 
equation  solving  technique  for  equation  29.  One  approach,  preferred  by 
some  in  the  finite  difference  community,  is  to  linearize  the  terms  in  the 
left  matrix  of  equation  29,  and  achieve  a  direct  solution.  The  alterna¬ 
tive  approach  is  to  employ  a  matrix  iteration  procedure  with  full  non¬ 
linearity  retained.  To  accomplish  the  former,  the  inverse  operation  [C]^ 
and  indicated  matrix  multiplications  in  equation  30  must  be  completed. 

Only  if  [C]  is  a  diagonal  matrix  is  the  inverse  trivial,  and  the  matrix  multi¬ 
plied  [u]k  and  [K]k  matrices  for  k  =  j,j+l,  become  modified  in  an  elementary 
manner,  say  [U]k  and  [K}k.  (Of  course,  this  diagonalizing  operation  destroys 
the  order-of-accuracy  intrinsic  to  the  basic  finite  element  algorithm  and  is 
to  be  avoided.)  Formally,  however,  the  solution  algorithm  for  equation 
29  can  be  made  non-iterative,  assuming  a  Taylor  series  expansion  for  the 

modified  matrix  terms  as,  e.g. 


I 

i 
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+  . . . 


(30) 


*  I®  K*  -  w,j' 

5[iii 

The  second  expression  is  obtained  assuming  the  dominant  non-linearity  has 
already  been  extracted  in  equation  26.  Then,  assuming  e  =  1/2  for 
simplicity,  equation  30  becomes  expressible  as 


which  is  a  linear  algebraic  equation  system  for  {Q> ^  and  directly  solvable. 

Based  upon  work  reported  herein,  the  steps  necessary  to  achieve  equation 
(31)  are  undesirable  in  terms  of  accuracy.  In  addition,  formation  of 
[C]k  and  [C]'1  are  to  be  avoided  at  all  cost.  While  [C]^  can  be  cleared 
from  the  left  of  equation  29  by  a  premultiplication  by  [C] j+ ^ ,  whereupon 
[C]j+1[l]  =  [C]j+1,  difficulty  appears  on  the  right  with  terms  in  [C]j+1[C]T’  / 
[i] .  This  is  removable  however,  by  limiting  0  =  1/2  which  produces  the 
trapezoidal  rule.  Formally,  using  equation  26  assembled  and  evaluated  at 
the  mid-interval,  equation  28  then  becomes  reexpressed  for  e  =  1/2  as 


No  difficulty  is  then  encountered  in  clearing  the  inverse,  and  equation  32 
takes  the  form 


By  assumption, 

* l  H*  *  t'lj)  <»> 
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and  the  Taylor  series  expansion  yields 


The  second  term  vanishes  identically  for  the  Euler  equations  in  conserva¬ 
tive  form,  and  the  linearized  direct  solution  form  for  equation  33  is 


Equation  36  retains  the  accuracy  of  the  basic  algorithm  and  is  uncondi¬ 
tionally  stable  according  to  a  linearized  analysis. 

The  desirable  feature  of  equation  36  is  that  it  is  linear,  hence 
solvable  in  a  non-iterative  manner  using  any  direct  procedure,  e.g.,  Gaussian 
elimination,  Cholesky  decomposition,  etc.  The  price  extracted  for  this 
feature  is  that,  while  the  [C] ,  [U] ,  and  [K]  matrices  are  banded,  they  be¬ 
come  block-banded  for  solution  of  systems  of  partial  differential  equations. 
Of  course,  direct  solutions  are  rarely  used  for  the  sake  of  computer  cost  and 
storage,  but  are  replaced  by  partial  factorization  of  the  matrix  (ref.  4). 

At  the  expense  of  extraneous  terms  being  added  to  the  differential  equation, 
this  operation  produces  block  tridiagonal  matrices  which  can  be  efficiently 
solved.  Hence,  the  order  of  accuracy  of  the  basic  algorithm,  equation  29, 
is  retained,  while  there  may  be  a  specific  degradation  of  accuracy,  level  at.  • 
any  given  solution  point. 

The  alternative  approach  is  to  matrix  solve  equation  33  as  non-linear. 
Letting  (F)  =  {0}  be  the  homogeneous  form  of  equation  33,  the  Newton 
iteration  algorithm  is 


The  dependent  variable  in  equation  37  is  the  iteration  vector,  defined 
by  the  relation 


+  6Q 


p+1 

j+1 


where  p  is  the  iteration  index.  The  right  side  of  equation  37  is 
equation  33  evaluated  with  the  p—  iterate,  i.e.. 


(F,Sn  ■  Se 


[c]e({Q)?tl  -  <q>j)  ♦  ^ge)?+ 


tge)j 


where 


(38) 


(39) 


(9e>S  =  (Me  +  [K]eHQ}£  +  (40) 


The  vanishing  of  { F}  to  within  definition  of  computed  zero  renders  equation 
37  homogeneous,  hence  convergence  of  the  iteration.  By  definition,  the 
Jacobian  is  the  derivative  of  equation  39  with  respect  to  Hence, 


M  -  sf 


[C]e+£  Me+  M, 


.  h 
+  7 


Mp  +  D<1 


3{Q}. 


(41) 


where  the  final  term  accounts  for  contributions  stemming  from  implicit  non¬ 
linearity.  The  rank  of  [j]  equals  the  order  of  {6Q};  Dirichlet  boundary 
constraints  are  applied  within  the  evaluation  of  {F}. 

The  ostensible  key  feature  of  the  iterative  solution  algorithm,  equa¬ 
tions  37  -  41  is  accuracy  and  elimination  (potentially)  of  the  block- 
banded"  Jacobian  matrix  structure  for  multiple  variable  systems.  Several' 
other  potentially  attractive  attributes  emerqe  on  closer  analysis.  First, 
as  discussed  for  the  non-iterative  algorithm,  three-dimensional  solution 
requirements  demand  a  spatial  factorization  that  permits  replacement  of  the 
large,  sparse-matrix  Jacobian  operations  with  elementary  banded-matrix 
procedures.  This  can  be  accomplished  in  the  present  instance  by  replacing 
the  multi -dimensional  interpolation  basis  {N^ ( 5t) } ,  equation  20,  with  the 
tensor  product  basis  (ref.  13) 

{NkU)>  =MNk(xa)}®  (Nk(x&)}®{Nk(xY)}  (42) 
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where  (x)  signifies  the  tensor  product.  The  Greek  letter  subscripts  are 
not  summation  indices  but  refer  instead  to  scalar  components  of  the 
resolution  of  x  on  to  the  coordinate  system  spanning  Rn.  Using  equation  42, 
the  original  finite  element  algorithm  statement,  equation  26,  becomes 
replaced  by  the  tensor  matrix  product  form.  For  example,  the  first  matrix 
becomes  of  the  form 

Me  ♦•Pde®  Me®  tc*le  M3) 

Hence,  the  Jacobian  can  be  evaluated  as  the  tensor  matrix  product, 

[J]  *=»  [Ji]@  [02](x)  [J3]  (44) 

and  each  of  the  [J  1  exhibit  the  desired  narrow  band  matrix  structure.  Spe- 

a 

cifically,  for  k  =  1  in  equation  42,  [J^]  is  tridiagonal,  while  for  k  =  2 
it  is  pentadiagonal  with  a  cyclic  profile. 

This  factorization  yields  a  highly  efficient  solution  algorithm  upon 
application  of  the  concept  of  fractional  steps  (ref.  14).  Specifically,  eval¬ 
uation  of  the  homogeneous  solution  statement  { F } ,  equation  39,  is  replaced 
by  the  tensor  product  form,  and  the  Newton  algorithm  becomes 

[MQ>Stl]  ®  [j>(Q)J+1]  ®  * 

-  [f,(Q)5+1J  (45) 

•  •*•••  *  . 

Here,  Q  and  (j  represent  intermediate  iterates,  and  <5Q  is  interpreted  as  the 
iteration  vector  for  each  respective  iterate.  If  required  for  accuracy, 
evaluation  of  the  right-side  of  equation  (45)  may  be  replaced  with  conven¬ 
tional  multi-dimensional  operations. 

4.  Accuracy  and  Convergence,  Stability 

A  von  Neumann  stability  analysis  (ref.  2)  can  quantify  the  formal  order 
of  accuracy  of  the  developed  algorithm  for  a  linearized  one-dimensional 

15- 

■  ,.i  j  m  m  w  +  ,i,  ■*>  *' 


I 


t 
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equation,  hence  predict  an  appropriate  value  for  £  equation  25.  There¬ 
fore,  consider  the  inviscid  Xj  momentum  equation  2  with  constant  ad- 
vection  velocity  U0< 


l»=i+Uo£*°  (46) 

The  analytical  solution  for  equation  46  is  the  Fourier  expansion 

u(x,t)  =  V  exp  Jj  gj(x  -  U0t)]  (47) 

where  i  =  /T,  a)  =  2tt/X  is  the  wave  number  where  X  is  wavelength,  and  V  is 
the  initial  velocity  distribution. 

A  rearrangement  of  equation  25  yields  the  desired  form  for  analysis. 
Using  a  Green-Gauss  theorem,  the  second  term  in  equation  25  can  be  trans¬ 
formed  to 


$*f  fNk>VL(qe)  =  -V.JtJ  {Nk}L(qej]  +  (Nk)L(qe)n 

R"  R"  8R _ 


(48) 


Since  the  tensor  product  algorithm  form  is  to  be  utilized,  from  reference  15 

let  $  h  v  ,  where  v  is  a  scalar  to  be  determined  and  a“  is  the  measure 
a  e  a  a  a  e 

of  the  discretization  of  the  one-dimensional  domain  Rg  spanned  by  xq  with 
unit  vector  ?  .  Under  this  assumption,  the  second  term  in  equation  48 
vanishes  identically.  Setting  X  =  0  for  convenience  for  this  analysis, 
equation  25  can  be  written  in  the  compact  tensor  product  form 


S 


e 


_  a  ‘ 

L(qe) 

{0} 


(49) 


It  is  an  elementary  exercise  to  evaluate  equation  49  for  any  differen¬ 
tial  equation.  Recalling  equation  26  as  the  required  final  form,  the 
initial-value  and  convection  matrices  are  defined  as 


(Nk(xa))Tdx 


a 


(50) 
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(51) 


Ptfe  2  Lit  +  5r[lU“^{Nk,'Nk,T]dX“ 

*e  L  J  L 

For  equation  46,  and  limitinq  the  interpolation  basis  to  linear,  i.e., 
k  =  1  in  equation  42,  equation  50  yields 


A“  [2  1 
^Je  =  5“  j  2 


«  Ot 

+  vAe 
+  T~ 


-1  -1 


for  v  assumed  a  constant.  For  a  uniform  discretization,  the  assembly  of 
this  term  in  equation  26  over  the  two  elements  sharing  node  j  yields  the 
finite  difference  recursion  relation 


se  "  f  (H3v)  Q'.t  *  4Q:  *  (l-3v)Q'+1  (53 

where  Aa  is  the  measure  of  the  uniform  discretization  of  Ra.  Similarly,  the 
convection  operator  in  equation  (46)  produces  (ref.  7) 


For  a  constant  advection  velocity,  {Ug>e  =  UQ { 1 } ;  the  assembly  yields  the 
finite  difference  recursion  relation 

Se  [UjeWg  -  ^-(l+ZvjQj.j  +  4vQ.  +  (l-2v)Qj+1 

=  T~  _Qj-l  +  Qj+1  +  V  _Qj-l  +  2Qj  ‘  Qj+1  (55) 

In  equation  55,  the  first  term  is  the  familiar  second  order  accurate 
central  difference  convection  operation.  The  second  term  is  equally  familiar 
as  the  second-order  accurate  diffusion  operator  with  UQvAa  as  the  "viscosity" 
level . 
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The  occurrence  of  these  two  difference  expressions  is  to  be  expected, 
as  confirmation  of  the  assertion  that  linear  finite  element  methodology 
yields  optimal  order-accuracy  difference  relations  for  the  same  degree 
interpolation.  The  full  impact  of  the  algorithmic  development  becomes 
apparent  upon  completion  of  the  von  Newman  stability  analysis.  For  the 
complete  algorithm  applied  to  equation  46,  formal  order-of-accuracy  is 
assessed  using  the  semi -discrete  Fourier  expansion 

u*(Ax,t)  e  zgue  =  V  exp  jj  a)(jAx  -  Xt)J  (56) 

where  Ax  is  the  measure  of  the  uniform  discretization  of  R1  and  X  e  3  +  i 6 , 
where  6  and  6  are  real  numbers.  Direct  substitution  of  equation  56  into 
equations  53  and  55,  and  expanding  the  resultant  expressions  for  p  and  <5 
in  a  Taylor  series  yields  (ref.  15) 


Here,  d  =  u>Ax  and  0(  )  indicates  the  order  of  the  truncated  term.  The  phase 
accuracy  of  the  discrete  solution  u*  agrees  with  UQ  to  sixth  order  in  Ax  by 
requiring  v"1  e  /IS".  For  v  >  0,  correspondingly  6  <  0,  and  an  artificial 
damping  is  introduced  of  the  form, 

exp jj-  to4kt  +  0(Ax5)J  (59) 

where  k  =  U0(Ax)3/12/IT  is  the  damping  coefficient.  The  damping  is  highly 
selective  due  to  the  w4  factor. 

Setting  v  e  0  eliminates  6,  and  the  derived  algorithm  is  a  fourth- 
order  accurate  representation  of  the  differential  equation  46.  This  high 
order  accuracy,  obtained  with  use  of  the  simplest  linear  interpolation,  re¬ 
sults  directly  from  the  finite  element-derived  initial-value  matrix  [Ca]e- 
In  distinction,  finite  difference  practice  replaces  equation  53  with  AaQ.' 

J 

which  immediately  degrades  the  algorithm  to  second  order  accuracy.  In 


addition,  the  conventional  finite  difference  approach  to  introduce 
dissipation  is  to  add  an  "artificial  viscosity"  term  p  92u/9x2  to 
equation  46.  The  resulting  dissipation  term  directly  comparable  to 
equation  59  is 

exp|~-w2pt  +  0(Ax2)J  (60) 

The  resultant  artificial  damping  is  less  selective  due  to  the  wave  number 
exponent  of  2  and  the  appearance  of  the  term  of  order  Ax2. 

Using  a  fully  discrete  Fourier  analysis,  see  reference  16,  the  ampli¬ 
fication  factor  for  the  dissipative  finite  element  algorithm  for  9  e  ^  is 


1+h  cos(<dAx)  -  3Cv  sin2(4aiAx)  -  i|-(C+2v)sin(u)Ax) 

■ - - ^ - 

1+H  cos(aiAx)  +  3Cv  sin2(%i)Ax)  +  i^{C-2v)sin(uAx) 


where  C  =  UQAt/Ax  is  the  Courant  Number.  The  numerator  and  denominator  are 
complex  conjugates  for  v  =  0,  hence  the  basic  (non-dissipative)  algorithm  is 
neutrally  stable  for  all  Ax  and  At.  Therefore,  error  induced  by  a  solution 
will  propagate  undamped  and  unmagnified  throughout  the  solution  domain. 
Selecting  v  >  0  destroys  the  neutral  stability,  and  damping  of  both  induced 
error  and  the  physical  solution  will  result. 


SECTION  III 


DISCUSSION  OF  RESULTS 


1.  Overview 

The  convergence  and  stability  analyses  are  strictly  limited  to  the 
elementary  differential  equation  studied.  In  particular,  the  practical 
instance  of  a  constant  imposed  velocity  is  negligible.  Nevertheless, 
they  do  help  interpret  the  predictions  of  numerical  solutions  regarding 
indicated  error  sources  and  their  modification.  Extension  of  the  elemen¬ 
tary  analysis  to  include  use  of  the  quadratic  basis  within  the  algorithm 
is  complicated  by  the  fact  that  a  simple  finite  difference  nodal  recursion 
relation  does  not  result.  Therefore,  the  analysis  of  carefully  executed 
numerical  tests  appears  the  appropriate  method  to  refine  and  extend  the 
scope  of  these  predictions.  For  this  purpose,  linear  and  non-linear  dif¬ 
ferential  equations,  that  model  the  convection  operator  intrinsic  to  the 
Navier-Stokes  equations,  have  been  selected  for  study.  Primary  emphasis  is 
an  assessment  of  factors  affecting  solution  accuracy  and  economy.  Key 
results  of  the  study  are  discussed  in  this  section. 

2.  Accuracy  and  Error  Control 

The  developed  finite  element  solution  algorithm  requires  evaluation 
with  respect  to  parameters  to  be  selected  for  error  control.  Specifically, 
these  include 

a)  discretization/integration  step  (Courant  Number),  C 

b)  finite  element  interpolation  degree,  k 

c)  phase  error  control,  v 

Recall  the  definition  of  Courant  Number,  C  e  UAt/Ag,  which  is  available  of 
course  in  any  numerical  algorithm.  The  last  two  are  more  uniquely  expressed 
within  the  structure  of  the  developed  algorithm. 

A  demanding  test  problem  corresponds  to  the  two  dimensional  form  of 
equation  46, 


where  =  rwB  corresponds  to  a  solid  body  rotation  with  angular  velocity 
a),  and  the  initial  distribution  for  u  is  the  "cosine  hill."  The  initial 
distribution  would  appear  identical  to  Figure  la),  moved  to  the  9  o'clock 
position,  and  the  exact  solution  is  diffusion-  and  dispersion-free  con¬ 
vection  clockwise  around  the  solution  domain.  From  reference  16,  Figures 
2a)  -  c)  illustrate  a  typical  linear  element  solution  at  the  quarter,  three- 
quarter  and  full -rotation  time  step,  for  the  non-dissipative  form  of  the 
algorithm.  The  peak  value  remains  within  1%  of  its  original  level  of  100, 
although  the  lagging  phase  error  retards  its  relative  displacement.  The 
ripple  structure  in  the  background  plane  is  pure  dispersion  error,  which 
propagates  generally  undissipated  throughout  the  solution,  in  agreement  with 
the  basic  neutral  stability.  The  results  shown  in  Figure  2d)  were  obtained 
using  a  digital  filtering  technique  that  effectively  dissipates  the  back¬ 
ground  error  structure.  However,  the  concentration  peak  was  also  diffused 
well  below  its  initial-value.  Hence,  these  solutions  are  quite  good  in 
comparison  to  standard  finite  difference  results,  yet  exhibit  significant 
error  upon  a  closer  examination. 

A  comprehensive  numerical  evaluation  of  this  problem,  for  variation  of 
parameters  a)  -c),  can  firmly  quantize  algorithm  performance  and  has  been 
completed.  For  comparison.  Figure  3  shows  the  appropriate  portion  of  the 
actual  initial  distribution  of  u,  which  also  corresponds  to  the  exact  final 
solution.  Figure  3b)  shows  the  distribution  computed  for  k  *  1,  v  =  0  and 
C  =  0.25  (at  the  center  of  the  solution  distribution;  since  UQ  varies 
linearly  with  r,  0  <  C  <  0.4  on  the  field).  As  noted,  the  peak  lags  by 
about  ^Ae  of  occurring  at  the  correct  node  point.  The  background  ripple 
structure,  observed  in  Figure  2,  occurs  as  nominal  rows  of  positive  and  nega- 
tive  concentrations  in  Figure  3b).  This  phase  error  induced  extremum  level 
is  -10  and  the  nominal  wavelength  is  4Ae>  As  a  consequence,  the  symmetries 
of  the  correct  solution  are  not  apparent.  However,  if  two  additional  time 
steps  were  taken,  the  symmetries  would  nominally  reappear  and  a  peak  value 
of  99  would  occur  at  the  circled  location.  Hence,  the  solution  visually 
appears  a  close  approximation  to  the  correct  solution.  Figure  2c). 

Increasing  the  integration  time  step  (or,  equivalently,  decreasing  the 
mesh  measure)  increases  the  Courant  number  which  adds  truncation  error  to 
the  numerical  solution.  Comparing  Figure  4b),  obtained  using  C  =  0.5,  k  =  1 
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Figure  2.  Advection  of  Cosine  Hill  in  a  Solid-Body 

Rotation  Velocity  Field,  C  =  0.1  (Reference  16) 
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Figure  3.  Advection  of  Cosine  Hill  in  Solid  Body  Rotation, 
Initial  Condition,  Exact  Final  Solution 
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a)  C  =  0.25,  v  =  0 
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b)  C  =  0.5,  v  =  0 


Figure  4.  Advection  of  Cosine  Hill  In  Solid  Body  Rotation, 
Linear  Tensor  Product  Algorithm 
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and  v  =  0,  to  Figure  4a)  shows  the  noticeable  degradation.  The  computed 
peak  lags  by  one  full  mesh  interval,  and  the  associated  dispersion  error 
has  reduced  its  level  to  93.  However,  the  solution  symmetry  about  the 
indicated  peak  is  quite  good,  (note  dashed  contour),  and  it  would  appear 
acceptable  visually.  The  trailing  dispersion  error  wake  extremum  is 
increased  to  -17,  and  the  wavelength  is  nominally  centered  on  5Ag.  For 
comparison.  Figure  4c)  shows  the  final  solution  as  predicted  using  the  Crank- 
Nicolson  finite  difference  algorithm  at  C  =  0.25.  The  result  is  much  poorer 
than  the  linear  element  tensor  product  solution  at  twice  the  Courant  number, 
i.e.,  half  the  computer  CPU.  The  peak  is  only  47  and  the  extremum  dispersion 
error  level  is  -24.  This  popular  algorithm  is  obtained  by  the  elementary 
change  of  equation  52  to 


3  O' 
0  3 


(63) 


As  seen,  the  consequence  is  introduction  of  unacceptable  error  levels  at  mod¬ 
est  Courant  Number. 

A  prime  requirement  is  to  evaluate  the  efficiency  of  v,  the  introduced 
dissipation  parameter,  to  alter  (improve)  the  error  levels  in  these  solutions. 
Elementary  one-dimensional  studies,  references  16  -  17,  have  confirmed 
that  the  optimal  linear  analysis  value  of  v  =  1//I3"  =  0.258  is  entirely  too 
large.  The  loop  test  case  is  much  more  demanding,  and  confirmed  the  viability 
of  separate  and  distinct  value  of  v,  for  use  in  the  initial-value  and  convec¬ 
tion  matrices,  see  equations  52  and  53.  Terming  them  correspondingly 
v1  and  vd,  Table  1  summarizes  the  results  of  the  cosine  hill  tests  for  the 
linear  tensor  product  algorithm.  All  data  are  for  the  final  station.  For 
comparison,  case  1  is  the  summary  for  the  solution  shown  in  Figure  4b, 
in  terms  of  computed  solution  peak,  displacement  (in  Afi)  of  the  peak  from  its 
correct  location,  and  fore/aft  levels,  hence  symmetry  about  the  peak.  The 
dispersion  wake  error  is  quantized  on  extremum  level  and  nominal  wavelength 
of  the  dominant  trailing  wave.  Using  vd  =  0.06  alone,  case  2,  which 
amounts  solely  to  an  artificial  diffusion,  the  wake  error  is  decreased  by 
two  while  the  solution  peak  is  decreased  by  15%  to  78.  Introducing  v"*  <  vd 
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Figure  4.  Advectlon  of  Cosine  Hill  In  Solid  Body  Rotation, 
Linear  Tensor  Product  Algorithm  (Concluded) 
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TABLE  1 


SUMMARY  OF  COSINE  HILL  SENSITIVITY  STUDY 
LINEAR  TENSOR  PRODUCT  ALGORITHM,  C  =  0.5 


Case 

Dissipation 

Solution 

Wake  Dispersion  Error 

No. 

v1 

vd 

Peak 

Displacement* ** 

Symmetry 

Peak 

Wavelength* 

1 

.0 

.0 

93 

-1 

81/82 

-17 

4 

2 

.0 

.06 

78 

-1 

71/66 

-9 

5 

3 

.012 

.06 

102 

-1 

87/89 

-21 

6 

4 

.012 

.072 

98 

-1 

85/86 

-19 

6 

5 

.015 

.06 

110 

-1 

92/98 

-26 

6 

6 

.012 

.09 

93 

-1 

81/80 

-16 

6 

7 

.03 

.06 

117 

-1 

61/161 

-74 

*Units  are  Ag,  the  uniform  discretization  measure. 

**No  discernable  dominant  wavelength. 
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over  a  range,  case  Numbers  3-7  can  remarkably  increase  the  computed 
solution  peak  and  affect  solution  symmetry  as  well  as  dispersion  error. 

Case  4  corresponds  to  the  "numerically  optimized"  combination,  and 
the  complete  final  solution  station  is  shown  in  Figure  4d  .  Comparing  to 
Figure  4b,  and  to  the  correct  solution  in  Figure  3,  the  v-  modified 
solution  is  a  uniform  improvement  except  for  the  modest  increase  in  wake 
error  level.  Importantly,  the  peak  of  98  and  the  surrounding  level  near  85 
are  both  substantial  improvements.  The  ability  of  a  small  proportion  of 
v1  to  totally  alter  the  solution,  and  the  pure  diffusive  character  of  vd, 
should  be  of  considerable  use  and  certainly  requires  a  definitive  study. 

The  same  test  sequence  has  been  completed  for  the  quadratic  tensor 
product  algorithm.  Figure  5a  is  the  base  solution  for  C  =  0.25  and  v  =  0. 
Comparing  to  Figure  3,  it  is  an  excellent  approximation  to  the  correct 
solution.  The  peak  occurs  exactly  at  the  correct  node,  and  the  surrounding 
solution  symmetry  and  levels  are  excellent.  Furthermore,  the  extremum  wake 
error  level  is  only  -2.  On  all  assessment  bases,  the  improvement  over  the 
linear  tensor  product  results,  see  Figure  4a,  and  for  the  Crank-Nicolson 
results  Figure  4c,  is  impressive.  Equally  impressive  is  the  observation 
that  only  1%  additional  computer  storage  was  required  and  the  increase  in 
CPU  was  only  16%  (due  mainly  to  a  pentadiagonal  rather  than  tridiagonal 
Jacobian).  Doubling  the  Courant  number  produced  the  results  shown  in  Figure 
5b),  which  are  a  nominal  improvement  over  the  linear  tensor  product  results 
at  C  =  0.25,  Figure  4i.  Hence,  comparable  accuracy  is  attained  at  a  40% 
decrease  in  CPU.  Finally,  Figure  5c  displays  the  results  obtained  using 
the  quadratic  diagonalized  initial-value  matrix  equivalent  of  Crank-Nicolson. 
The  results  are  uniformly  poorer  on  all  bases  of  comparison. 

Table  2  summarizes  the  results  of  the  v  sensitivity  study  for  the  quad¬ 
ratic  tensor  product  algorithm.  Since  the  non-dissipative  results  maintained 
a  larger  peak  value,  the  algorithm  can  withstand  somewhat  larger  levels  of 
vd  in  comparison  to  the  linear  element  form.  Case  5  results  were  the 
"best"  improvement  and  the  entire  solution  is  shown  in  Figure  5d).  In 
comparison  to  the  linear  element  case,  the  improvement  is  not  quite  as 
substantial  (since  the  error  requiring  improvement  was  smaller),  but  never¬ 
theless  uniformly  positive. 
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Figure  5.  Advectlon  of  Cosine  Hill  In  Solid  Body  Rotation, 
Quadratic  Tensor  Product  Algorithm 
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d)  C  =  0.5,  v1  =  0.012,  vd  =  0.09 


Figure  5.  Advection  of  Cosine  Hill  in  Solid  Body  Rotation, 
Quadratic  Tensor  Product  Algorithm  (Concluded) 
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TABLE  2 

SUMMARY  OF  COSINE  HILL  SENSITIVITY  STUDY 
QUADRATIC  TENSOR  PRODUCT  ALGORITHM,  C  =  0.5 


Comparison  sensitivity  studies  have  been  essentially  completed  for 
the  non-linear  form  of  the  model  equation  62,  i.e., 

L(u)  =  |£  +  T5*Vu  =  0  (64) 

where  tJ  =  ui  +  vj  and  v  satisfies  an  identical  equation.  The  test  case  is 
the  two-dimensional  square  wave  propagating  parallel  to  the  principal 
diagonal  of  the  solution  domain.  Figure  6a  shows  the  initial  condition, 
while  Figure  6b  shows  the  final  solution  obtained  for  k  =  1,  v  =  v1  = 

Vd  =  0.182  and  C  =  0.125.  This  solution  is  a  close  approximation  to  the 
exact  solution,  wherein  the  dominant  wave  front  remains  interpolated 
across  only  one  element,  the  trailing  plateau  remains  at  unity  with  no 
overshoot,  and  the  solution  gradient  perpendicular  to  the  left  boundary 
vanishes.  Increasing  vd  >  0.182  levels  the  indicated  overshoots  but  also 
diffuses  the  wave  front  across  three  or  more  A  .  The  action  of  a  separately 
evaluated  v1  ^  vQ  in  correcting  this  trend  needs  to  be  evaluated.  Figure  6c 
is  the  comparison  finite  difference  solution,  which  is  essentially  identical 
to  that  in  reference  4,  page  96.  It  is  obtained  by  use  of  equation  63  instead  of 
equation  52,  and  replacement  of  the  v-modified  term  in  equation  54  with 
an  artificial  viscosity  term.  The  influence  of  the  degraded  selectivity  of 
the  diffusion  operator,  compare  equations  59  and  60,  is  quite  apparent. 

Again,  the  interpolation  theory-derived  algorithm  appears  superior  to  the 
equivalent  complexity  finite  difference  algorithm. 

The  action  of  v  >  0  within  the  dissipative  Weighted-Residuals  theory 
is  crucial.  Specifically,  Figure  7a  shows  the  solution  generated  by  the 
non-dissipati ve  k  =  1  algorithm,  which  has  become  completely  nonsensical. 

Hence,  it  is  imperative  that  the  basic  neutral  stability  be  modified  by 
v  >  0  for  a  non-linear  problem  involving  steep  field  gradients,  e.g.  ,  shocks. 

The  same  holds  for  the  quadratic  embodiment  of  the  algorithm.  Figure  7b 
shows  the  final  solution  obtained  for  k  =  2,  vd  =  0.182,  =  0  and  C  =  0.125. 

Solution  fidelity  is  comparable  to  that  of  the  k  =  1  solution,  with  a  slight 
improvement  in  steepness  of  the  dominant  wave  front.  A  sensitivity  study  for 
v*  >  0  should  also  be  completed. 

F 


a)  Initial  Condition,  t  =  0 


b)  nAt  =  120,  v1  =  vd  =  0.182 


c}  nAt  =  120,  Finite  Difference  Form 


Figure  6,  Two-Dimensional  Square  Wave  Convection, 

Linear  Tensor  Product  Algorithm,  C  =  0.125 


3.  Matrix  Iterative  Solution  Efficiency 

The  solutions  discussed  were  generated  employing  the  Newton  iteration 
algorithm,  equations  37  -  41,  as  embodied  within  the  tensor  product 
formulation,  equation  45.  This  matrix  solution  algorithm  can  be  rendered 
non-iterative  by  the  elementary  modifications  to  (F> ,  equation  39,  of 
changing  the  sign  of  h/2  and  removing  all  terms  involving  matrix  products  on 
{Q}j+1-  The  Jacobian  is  unaltered.  Most  of  the  linear  cosine-hill  convec¬ 
tion  problems  were  solved  using  the  non-iterative  form.  With  a  con¬ 
vergence  requirement  of  10’ 3 ,  the  iterative  algorithm  required  two  iterations 
to  produce  solutions  essentially  identical  to  the  direct  solution.  Hence,  a 
CPU  savings  of  50%  can  be  expected  for  a  single  scalar  equation.  The  intro¬ 
duction  of  multiple  dependent  variables  would  introduce  a  block-banded 
structure  into  [J]  that  would  moderate  this  basic  advantage.  The  iterative 
form  would  employ  multiple  right-side  back-substitutions  into  the  LU  de¬ 
composition  of  the  basic  tri-  or  pentadiagonal  Jacobian  [j] .  At  some  point 
a  trade-off  might  occur,  and  the  test  problem  must  correspond  to  solution 
of  multiple  equations. 

In  anticipation,  the  fully  iterative  non-linear  algorithm  has  been 
evaluated  for  the  turbulent  boundary  layer  equation  using  both  mixing 
length  and  turbulence  kinetic  energy-dissipation  function  turbulence  clo¬ 
sures.  The  fundamental  question  was  whether  an  economical  prediction  tech¬ 
nique,  based  upon  finite  differences,  could  yield  a  reduced  number  of 
iterations  for  this  rather  non-linear  equation  system.  One  procedure  was 
the  non-uniform-step  extension  of  the  finite  difference  method  used  by 
Soliman  (ref.  17).  A  second  order  accurate  estimate  of  qj+j  is  obtained 
using  the  leap-frog  finite  difference  equation  for  non-uniform  grid 

qj+l  =  qj  t1  ’  +  n)  ♦  (1  +  n)z]  +  [1  +  nl^j+i 

+  9j_i  [?(!  +  n)  -  (1  +  n)2J  (65) 

where  n  =  A is  the  non-uniformity  of  the  integration  step.  A  second- 
order  accurate  difference  formula  for  q^  =  qj(qj,  9j_i>  9j_2^  was  employed. 
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Hence,  this  approach  constitutes  a  second-order  accurate  leap-frog  pre¬ 
dictor,  followed  by  the  second-order  accurate  trapezoidal  integration 
corrector.  The  intent  is  to  accurately  predict  q^+1  such  that  fewer 
corrections  are  needed  for  convergence.  No  equation  solution  is  required 
but  computer  CPU  cost  is  incurred  from  the  finite  difference  evaluations. 


The  second  predictor  technique  is  conceptually  opposite,  in  that  no 
attempt  is  made  to  predict  qj+^  but  instead  a  maximum  evaluation  of  6qj+1 
is  sought.  Assuming  the  second  and  fourth  terms  in  equation  39  are 


stored  separately,  and  further  that  (f}j+j 

ql+1  =  q^  eliminates  the  first  term  in  equation  39  and  the  estimate  of 


{f}.,  then  specifying  that 

J 


{F}j+i  is  constructed  directly  on  q^.  This  method  is  termed  <5Q-predict, 
requires  a  matrix  solution  of  equation  37  but  eliminates  the  element  DO 
loop  to  form  (F)j+p  equation  39.  The  comparison  method  is  straight¬ 
forward  application  of  the  Newton  algorithm.  As  with  <$Q-predict,  qj+1  = 
but  a  computer  program  pass  is  made  to  correctly  evaluate  {g}1.  .,  hence 
<F(Q‘.n)>. 


The  results  are  summarized  in  Table  3  and  categorized  with  respect  to 
predictor  type  and  input  solution  parameters.  Decreasing  the  convergence 
requirement  to  10”1*,  Tests  1-2,  reduced  the  number  of  {F}  evaluations  and 
matrix  solves  by  about  40%  and  did  not  measurably  affect  accuracy,  as  measured 
in  the  shape  factor  norm.  The  accuracy  of  both  predictor-type  procedures 
was  identical  to  the  Newton  solution  as  measured  in  both  the  energy  and 
shape  factor  norms.  Test  2.  Surprisingly,  the  Q-predict  using  the  finite 
difference  formula  did  not  reduce  the  number  of  F  evaluations,  i.e.,  {Q> j +1 
could  not  be  predicted  with  sufficient  accuracy.  As  a  consequence,  its 
solution  CPU  cost  was  maximum.  The  <5Q-predict  required  fewer  solution  steps 
and  {F}  evaluations,  but  more  equation  solutions  (207  +  80  =  287  vs  272)  than 
the  Newton  algorithm  for  a  net  13%  savings.  Increasing  the  maximum 
T-step  by  a  factor  of  six  does  not  change  the  overall  solution  accuracy, 
compare  tests  2-3,  nor  the  overall  comparison  among  the  solution  methods. 
Hence,  the  Q-predict  method  can  be  discarded.  For  test  4,  the  maximum  T-step 
was  reached  earlier  in  the  solution,  and  the  increase  in  truncation  error 
has  decreased  the  solution  energy  by  about  a  percent.  Hence,  the  mathemati¬ 
cal  measure  of  the  solution  is  degraded,  but  the  engineering  measure  of 
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shape  factor  norm  is  negligibly  different.  Here,  the  6Q-predict  was 
optimized  to  allow  maximal  early  steps  which  reduced  the  total  solution 
steps  by  seven  to  30.  However,  the  increased  number  of  matrix  solutions 
(30  +  111  =  141)  reduced  solution  efficiency  improvement  to  only  5%. 

Enhancing  the  no-predict  procedure  could  reduce  this  modest  superiority 
to  zero;  hence  the  6Q-predict  method  is  only  marginally  attractive.  There¬ 
fore,  neither  prediction  method  holds  particular  promise  for  economizing  an 
iterative  algorithm,  and  the  storage  and  logic  required  for  their  execution 
can  be  deleted. 

4.  Resolution  of  Gradients,  Shocks 

An  important  requirement  is  the  capturing  of  imbedded  flowfield  shocks. 
The  performance  criteria  for  an  algorithm  include  shock  steepness,  strength, 
and  location.  By  and  large,  the  addition  of  artificial  diffusion  to  an 
Euler  algorithm  permits  prediction  of  interpolated  approximations  of  a  shock. 
Use  of  conservative  or  non-conservative  forms  of  an  algorithm  affects  shock 
strength. 

The  results  shown  in  Figures  6-7,  for  the  two-dimensional  Burger's 
equation,  indicate  differences  in  travelling  wave  front  steepness  for  two 
forms  of  the  tensor  product  algorithm  in  comparison  to  a  finite  difference 
algorithm.  Additional  insight  can  be  obtained  from  examination  of  a  solu¬ 
tion  of  the  Euler  equations  for  one-dimensional  shocked  duct  flow.  The 
axial  momentum  equation  in  non-conservative  form  is 

<66> 

where  A(x)  is  the  duct  cross-sectional  area  distribution.  Figure  8a  shows 
the  initial  distribution  of  u(x)  with  a  shock  located  in  a  region  of  fine 
discretization;  each  symbol  corresponds  to  a  node  point.  Figure  8b  shows 
the  u  distribution  resulting  from  application  of  a  perturbation  to  p(x)  that 
moves  the  shock  a  modest  distance  upstream. 

The  influence  of  the  levels  of  the  dissipation  parameters  v’  and  vd, 
on  the  solution  of  u(x,t)  in  the  immediate  vicinity  of  the  perturbation  have 
been  studied.  Particular  interest  focuses  on  the  dispersion  error  amplitude 


Figure  8.  Initial  and  Pressure-Perturbed  velocity  uisii 
One-Dimensional  Supersonic  Shocked  Duct  Flow 


on  the  supersonic  side  of  the  shock,  and  on  the  decay  of  the  downstream 
traveling  wave.  Figure  9a  shows  the  preponderance  of  error  generated  by  the 
non-dissipative  form  of  the  algorithm  after  twenty  time  steps.  In  contrast, 
for  v1  =  vd  =  0.0645,  Figures  9b  -  9c  show  the  linear  element  solutions  at 
20  and  200  time  steps.  The  dispersion  error  waveform  on  the  supersonic 
side  remains  essentially  stationary  throughout,  and  the  amplitude  of  the 
travelling  wave  decays  rapidly. 

The  final  shock  strength  is  overpredicted  by  about  10%,  and  depends 
primarily  on  the  level  of  vd.  Figure  10  illustrates  the  influence  of  various 
vd  for  v1  =  0  for  the  linear  element  algorithm  at  nAt  =  20.  In  comparison  to 
Figure  10a,  which  is  the  correct  solution,  and  the  dashed  curves,  decreasing 
vd  in  the  range  0.2582  >  vd  _>  0.0645  improves  shock  strength  accuracy, 
while  admitting  a  larger  dispersion  error  envelope  on  the  supersonic  side  and 
a  diminished  dissipation  of  the  downstream  travelling  wave.  Introduction  of 
v1  >  0  has  very  little  effect  on  shock  strength  or  dispersion  error  envelope, 
but  does  introduce  destructive  interference  into  the  travelling  wave.  Figure 
11  illustrates  the  results  obtained  for  the  linear  element  algorithm  for 
vd  =  0.0645,  0.0322  £  v1  £0.193.  Comparing  Figure  lid  with  lOd,  it  appears 
that  a  further  decrease  in  vd  to  improve  shock  strength  would  be  permissible. 

5.  Transformation  to  Arbitrary  Domains 

For  the  tensor  product  algorithm  to  be  viable,  it  must  be  operable  for 
general  geometries.  Geometric  versatility  was  the  original  key  feature  of 
the  finite  element  concept,  but  use  of  coordinate  transformation  to  map  non¬ 
regular  shaped  domains  onto  the  unit  square  will  supplant  this.  Thompson 
and  co-workers,  references  18  -  20,  pioneered  the  concept  of  solution  of 
non-self-adjoint  Poisson  equations  to  generate  the  mappings.  An  alternative 
procedure,  using  an  analytical  transformation  on  a  local  basis,  appears  most 
useful  for  the  tensor  product  algorithm,  while  enjoying  direct  extension  to 
three-dimensional  space. 

Figure  12  illustrates  the  concept  for  a  straight-sided  domain  on  two- 
dimensional  space.  The  origin  of  the  "natural"  n^  coordinate  system  is 
defined  at  the  centroid.  The  coordinate  curves  of  n^  remain  straight  lines, 
but  the  system  is  not  orthogonal  in  the  plane.  However,  the  n^  system 
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Computed  Duct  Velocity  Distributions  Using  Linear  Element  Algorithm 


Figure  11.  Effect  of  Dissipation  Level  v  on  Shocked  Duct  Flow  Solution,  k  =  1,  v  =  0.0645,  nAt 


is  defined  as  orthogonal  and  normalized  in  the  transform  plane.  The  tool 
for  accomplishing  this,  on  a  local  basis,  is  contained  within  the  standard 
finite  element  isoparametric  coordinate  transformation  {Ni+  (n^)},  see 
references  12  and  21,  where 

(1  -  m)(l  -  n 2) 

1  J ( 1  +  n 1 ) ( 1  -  n2)l 

(Ni  (n.-)}  =  t  1(1  +  ni)(l  +  m)f  (67) 

+  1  4  (1  -  m)(l  +  tu) 


Employing  the  classical  concept  of  finite  element  interpolation  theory,  the 
global  x.  coordinate  system  is  interpolated  on  R*  in  terms  of  the  global  node 
coordinates  (XI }g  of  the  finite  element  domain  as 

Xi  =  (N1+  (ni ))T{Xl>e 

x2  =  (Mi+  (ni)}T(X2}e  (68) 

The  tensor  product  algorithm  requires  establishment  of  the  vector 
gradient,  equations  1-3.  Hence,  using  the  chain  rule 

3frqe(xi)  =  ^i)}T{Q}e 

J  J 
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(71) 


M 


a ( 1  -  q2)  +  b(l  +  ri2 ) 
a(l  -  m)  +  8(1  +  m) 


C ( 1  -  ri2 )  +  d(l  +  ri2 ) 
y(1  -  ni)  +  <5(1  +  m) 


where  a-d  and  a-6  are  specified  constants  that  depend  solely  on  the  global 
node  coordinates 


Equation 


From  equation 

M'1=  dirpj 


69  requires  the  inverse  operation,  ie. 

3rtj  3n2 

3xt  3Xi 
3n1  3n2 

8X2  8X2 

71,  and  the  definition  of  a  matrix  inverse, 
y(1  -  m)  +  <5(1  +  m)  »  -a(l  -  m )  -  8(1  +  m) 
-c(l  -  02)  -  d(l  +  ru)  .  a(l  -  n2)  +  b(l  +  112) 


(72) 


(73) 


where  det  [j]  is  the  determinant  of  the  Jacobian,  equation  71.  Hence,  since 
a-d  and  a-6  are  known  constants,  equation  69  is  formally  evaluated  as 


3jK(N>  -  tNl*  <74> 

J  ^ 

The  finite  element  algorithm  requires  integral  operations,  see  equation 
25.  The  differential  element  becomes 

""A 

■  det  [J]  dnjdrij  (75) 


For  the  single  vector  operation  within  the  tensor  product  algorithm,  det  [J] 
in  equations  75  and  73  will  cancel,  leaving  [J]”1  in  the  form  of  poly 
nomials.  The  unit  vector  Ej  is  contracted  with  the  convection  velocity; 
the  inner  product  is  an  invariant  and  should  be  evaluated  in  the  local  r^. 
coordinate  system. 

No  fundamental  change  is  involved  in  going  to  curved-sided  or  three- 
dimensional  element  domains.  The  extension  In  concept  is  illustrated  in 
Figure  13  for  a  two-dimensional  general  domain  the  x^  coordinate  system  is 
interpolated  quadra ti cal ly  as 
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Figure  12. 


(76) 


*1  =  {n2+  <ni )}T(X1 )e 

*2  =  (N2*  (n1 )}TCX2)e 


where  {XI}  contains  eight  elements  equal  to  the  respective  global  nodal 
coordinates.  The  matrix  elements  of  the  isoparametric  quadrilateral  basis 
are  expressed  in  the  ni  coordinate  system  as  (ref.  21). 


(1  -  n j ) ( 1  -  n2)(-n! 

(1  +  nj)(1  -  n2)(  n! 

(1  +  Hi ) (1  +  n2)(  Hi 

(1  -  nj ) ( 1  +  n2)(-nl 

2(1  -  nf )(1  -  n2) 

2(1  +  ru )<1  -  n l) 

2(1  -  nf)(l  +  n2) 

2(1  -  n,)0  +  n2) 


-  ^2 
-  r'2 
+  n2 
+  n2 


1) 

1) 

1) 

1) 


Then 


MV  -  SjM' 


(n)}T{Q}e 


(77) 


(78) 


and  the  number  of  matrix  elements  in  {Q}g  now  equals  eight.  Numerical 
integration  procedures  will  undoubtedly  be  required  to  evaluate  equation 
25  in  the  tensor  product  form  using  equation  78. 
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SECTION  IV 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  research  project  was  to  develop  a  highly  accurate 
and  efficient  numerical  solution  algorithm  for  the  non-linear  three-dimen¬ 
sional  Navier-Stokes  equations  in  aerodynamics  applications.  A  candidate 
algorithm  has  been  derived  employing  finite  element  interpolation  theory, 
the  non-linear  extension  of  quasi-variational  principles,  and  the  concept 
of  tensor  product  bases.  The  resultant  solution  statement  is  rendered 
soluable  using  an  implicit  integration  algorithm,  and  the  replacement  of 
the  standard  Jacobian  of  a  Newton  iteration  algorithm  with  a  tensor  matrix 
product  form. 

A  linearized  stability  analysis  indicates  the  basic  algorithm  is  spa¬ 
tially  fourth-order  accurate  in  its  most  elementary  embodiment.  The  control 
of  an  added  dissipation  mechanism  can  elevate  this  to  sixth  order,  but 
numerical  experimentation  indicates  the  resultant  artificial  diffusion  is 
unacceptably  large.  By  the  same  measure,  this  additional  accuracy  is  in¬ 
trinsic  to  the  quadratic  element  embodiment  of  the  algorithm  with  freedom 
from  artificial  diffusion.  The  results  of  several  numerical  experiments 
for  single-  and  multi -dimensional  convection-dominated  test  problems  confirms 
the  basic  viability  of  the  developed  algorithm  and  its  tensor  product  formu¬ 
lation.  The  latter  is  of  paramount  importance  in  rendering  the  algorithm 
nominally  as  efficient  as  familiar  lower-order  accurate  finite  difference 
formulations. 

These  results  of  the  first  year's  research  on  the  topic  are  highly 
encouraging.  Specific  recommendations  include: 

1.  A  computer  program  test  bed  should  be  constructed  for  solution 
of  the  mul ti -dimensional ,  multi -dependent  variable  Navier-Stokes 
equations.  This  program  should  contain  both  the  linear  and 
quadratic  element  embodiments  for  solution  comparisons. 

2.  Evaluation  of  the  transformation  of  arbitrarily  shaped  solution 
domains  to  a  regular  grid,  using  the  isoparametric  mapping,  should 
be  accomplished  in  this  program. 

3.  The  developed  numerical  procedures  should  be  evaluated  for  solu¬ 
tion  of  several  Navier-Stokes  problems  for  which  accurate  com¬ 
parison  results  are  available. 
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